%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mapping reduced form coefficients into structural coefficients
% Identification of investment-specific,
% skill-biased and unskill-biased technology shocks using
% long-run sign restrictions: here on a trivariate submatrix
% this version: May 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inputs:
% BB                reduced form coefficients for ndraws
% Omega             VCV from reduced form residuals for ndraws
% nlags             number of lags
% ndraws            number of draws
% HORIZON_ident     identification horizon
% c                 choice of constant in the reduced form VAR
% num_vars          number of variables in the VAR
% nshocks           number of identified shocks
% hndx              index of hours in the VAR
% do_level          hours in levels of first differences

% output:
% Response          Responses to one standard deviation of all shocks for all horizons and all draws
% Atilde            matrix A from the structural representation
% Respres           Response to (all) one standard deviation residual shock(s)
%                   size is num_vars,HORIZON_ident,ndraws

function [Response,Respres]=ident_sign_trijoint(BB,Omega,nlags,ndraws,HORIZON_ident,c,num_vars,nshocks,hndx,do_level)

Respdraws = [];
Respresdraws = [];
Adraws = [];

maxdraws = 10; % number of draws for the rotations
display('how long does one rotation last with:');
display(maxdraws);
tic
for n = 1:ndraws
    
    % Iteration
    dnum=num2str(n);
    if dnum(end)==num2str(0) && dnum(end-1)==num2str(0)
        disp('iteration'),disp(n)
    end
    
    % Step I:
    % transform the reduced coefficients into structural coefficents Phi
    Btilde = zeros(num_vars,num_vars,HORIZON_ident);
    for k=1:nlags
        Btilde(:,:,k) = BB(:,((k-1)*num_vars+c+1):(k*num_vars+c),n);
    end;
    Phi = zeros(num_vars,num_vars,HORIZON_ident+1);
    Phi(:,:,1) = eye(num_vars,num_vars); % this is time zero, i.e. the impact response
    
    for s = 2 : HORIZON_ident,
        Phitemp = zeros(num_vars,num_vars);
        for j= 1:(s-1),
            Phitemp = Phitemp + Phi(:,:,(s-j))*Btilde(:,:,j);
        end;
        Phi(:,:,s) = Phitemp;
    end;
    
    % Step II: Implementing the restrictions
    % i.e. get CCtilde and A, the identification horizon is flexible here
    % and may be adjusted
    CCtilde = eye(num_vars,num_vars);
    for j= 2:HORIZON_ident+1
        CCtilde = CCtilde + Phi(:,:,j);
    end;
    
    % we are now looking for some matrix A, such that AA'=Omega and such that
    % CCtilde*A is lower triangular
    Sigma=CCtilde*Omega(:,:,n)*CCtilde';
    Ifactor = chol(Sigma)';
    
    % second, apply sign-run restrictions to subsystem (2:end,2:end)
    % here, we need to 'correct' the matrix Omega
    sub_ndx = num_vars-1;
    CorrMat = zeros(sub_ndx,sub_ndx);
    for i=1:sub_ndx
        CorrMat(i,:) = Ifactor(2:end,1)'.*Ifactor(i+1,1);
    end
    Sigma_sub = Sigma(2:end,2:end)-CorrMat;
    
    % We will be decomposing the upper left of the long-run variance
    Slong_sign = Sigma_sub(1:3,1:3);
    % starting value for the decomposition and rotation can be any matrix
    % that fulfills orthogonality assumption, usually Cholesky factor of the VCV is taken...
    T=chol(Slong_sign)';
        
    % grid
    gridspace = 2*pi/(maxdraws-1);
    grid=zeros(maxdraws,1);
    for k=1:(maxdraws-1)
        grid(k+1,1) = grid(k,1)+gridspace;
    end
    
    % starting the rotations
    for dr = 1:maxdraws
        
        % rotation matrix
        theta1 = grid(dr);
        for drr = 1:maxdraws
            theta2 = grid(drr);
            for drrr = 1:maxdraws
                sep = 0;
                theta3 = grid(drrr);
                Q1 = eye(3);
                Q2 = eye(3);
                Q3 = eye(3);
                Q1(1,1) = cos(theta1);
                Q1(1,2) = -sin(theta1);
                Q1(2,1) = sin(theta1);
                Q1(2,2) = cos(theta1);
                Q2(2,2) = cos(theta2);
                Q2(2,3) = -sin(theta2);
                Q2(3,2) = sin(theta2);
                Q2(3,3) = cos(theta2);
                Q3(1,1) = cos(theta3);
                Q3(1,3) = -sin(theta3);
                Q3(3,1) = sin(theta3);
                Q3(3,3) = cos(theta3);
                Q = Q1*Q2*Q3;
                
                % rotated matrix of long-run effects
                L_sign = T*Q;
                
                % check sign restrictions
                % first shock is the supply shock
                % check first column first
                supshock = 1;
                sbtshock = 2;
                ubtshock = 3;
                z1 = [L_sign(1,supshock)>0 L_sign(2,supshock)<0];
                zz1=all(z1); % true if all elements in z1 are nonzero
                if zz1==1
                    z2 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                    zz2=all(z2); % true if all elements in z1 are nonzero
                    if zz2==0
                        sbtshock = 3;
                        ubtshock = 2;
                        z3 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                        zz3=all(z3);
                    end
                elseif zz1==0
                    supshock = 2;
                    sbtshock = 1;
                    ubtshock = 3;
                    mz1 = [L_sign(1,supshock)>0 L_sign(2,supshock)<0];
                    mzz1=all(mz1);
                    if mzz1==1
                        mz2 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                        mzz2=all(mz2); % true if all elements in z1 are nonzero
                        if mzz2==0
                            sbtshock = 3;
                            ubtshock = 1;
                            mz3 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                            mzz3=all(mz3);
                        end
                    elseif mzz1==0
                        supshock = 3;
                        sbtshock = 1;
                        ubtshock = 2;
                        nz1 = [L_sign(1,supshock)>0 L_sign(2,supshock)<0];
                        nzz1=all(nz1);
                        if nzz1==1
                            nz2 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                            nzz2=all(nz2); % true if all elements in z1 are nonzero
                            if nzz2==0
                                sbtshock = 2;
                                ubtshock = 1;
                                nz3 = [L_sign(1,sbtshock)>0 L_sign(2,sbtshock)>0 L_sign(3,sbtshock)>0];
                                nzz3=all(nz3);
                            end
                        end
                    end
                end
                % check UBT shock last
                uz = [L_sign(1,ubtshock)>=0 L_sign(2,ubtshock)>=0 L_sign(3,ubtshock)<=0];
                uzz=all(uz); % true if all elements in z1 are nonzero
                if uzz == 1
                    if zz1 == 1
                        if zz2 == 1 || zz3 == 1
                            sep = 1;
                        end
                    elseif mzz1 == 1
                        if mzz2 == 1 || mzz3 == 1
                            sep = 1;
                        end
                    elseif nzz1 == 1
                        if nzz2 == 1 || nzz3 == 1
                            sep = 1;
                        end
                    end
                end
                    
                L_sign = T*Q;
                
                % calculating the remainder of L
                L_sub = zeros(num_vars-1,num_vars-1);
                L_sub(1:3,1:3) = L_sign;
                % calculating the remaining elements of the frist two columns
                for ii = 4:num_vars-1
                    L_sub(ii,1:3) = Sigma_sub(1:3,ii)'/L_sign';
                end
        
                % decomposing the remainder with Cholesky (need to correct the long-run
                % VCV)
                sub_ndx = num_vars-1-3;
                CorrMat_sub = zeros(sub_ndx,sub_ndx);
                for i=1:sub_ndx
                    CorrMat_sub(i,:) = L_sub(4:end,1)'.*L_sub(i+3,1) + L_sub(4:end,2)'.*L_sub(i+3,2) + L_sub(4:end,3)'.*L_sub(i+3,3);;
                end
                Slong_subsub = Sigma_sub(4:num_vars-1,4:num_vars-1)-CorrMat_sub;
                L_subsub = chol(Slong_subsub)';
                L_sub(4:end,4:end) = L_subsub;

                L=zeros(num_vars,num_vars);
                L(:,1) = Ifactor(:,1);
                L(2:end,2:end) = L_sub;

                % calculating Atilde
                A_sign = CCtilde\L;
                
                % save the good draws
                if sep == 1
                    
                    Avec = reshape(A_sign,1,num_vars*num_vars);
                    Adraws = [  Adraws
                        Avec ];
                    
                    % Step III: IRFs for all shocks (separately)
                    Resp = zeros(num_vars*num_vars,HORIZON_ident);
                    for shock_counter = 1:num_vars,
                        shock_vector = zeros(num_vars,1);
                        if shock_counter == 2
                            shock_vector(supshock+1) = 1;
                        elseif shock_counter == 3
                            shock_vector(sbtshock+1) = 1;
                        elseif shock_counter == 4
                            shock_vector(ubtshock+1) = 1;
                        else
                            shock_vector(shock_counter) = 1;
                        end
                        CCneutilde = zeros(num_vars,num_vars);
                        for time_counter = 1 : HORIZON_ident,
                            CCneutilde = CCneutilde + Phi(:,:,time_counter);
                            Resphelp1 = CCneutilde*A_sign*shock_vector; % for differences
                            Resphelp2 = Phi(:,:,time_counter)*A_sign*shock_vector; % for levels
                            Resp((shock_counter-1)*num_vars+1:(shock_counter-1)*num_vars+num_vars,time_counter) = Resphelp1;
                            if do_level==1
                                Resp((shock_counter-1)*num_vars+hndx,time_counter) = Resphelp2(hndx,:)*100;
                            end
                        end;
                    end;
                    
                    Respdraws = [  Respdraws
                        Resp ];
                    
                    % IRFs for all residual shock (bunched together)
                    Resresp = zeros(num_vars,HORIZON_ident);
                    shock_vector = zeros(num_vars,1);
                    shock_vector(nshocks+1:num_vars) = 1;
                    CCneutilde = zeros(num_vars,num_vars);
                    for time_counter = 1 : HORIZON_ident,
                        CCneutilde = CCneutilde + Phi(:,:,time_counter);
                        Resphelp1 = CCneutilde*A_sign*shock_vector; % for differences
                        Resphelp2 = Phi(:,:,time_counter)*A_sign*shock_vector; % for levels
                        Resresp(:,time_counter) = Resphelp1;
                        if do_level==1
                            Resresp(hndx,time_counter) = Resphelp2(hndx,:)*100;
                        end
                    end;
                    
                    Respresdraws = [  Respresdraws
                        Resresp ];
                    
                end
            end
        end
    end
end; % draws loop
toc

gdraws = size(Adraws,1);
Atilde = zeros(num_vars,num_vars,gdraws);
Rsize = num_vars*num_vars;
Response = zeros(Rsize,HORIZON_ident,gdraws);
Respres = zeros(num_vars,HORIZON_ident,gdraws);
for g=1:gdraws
    % reshaping back A
    Atilde(:,:,g) = reshape(Adraws(g,:),num_vars,num_vars);
    % reshaping back the responses
    Response(:,:,g) = Respdraws((g-1)*Rsize+1:(g-1)*Rsize+Rsize,:);
    Respres(:,:,g) = Respresdraws((g-1)*num_vars+1:(g-1)*num_vars+num_vars,:);
end